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I. INTRODUCTION 



The heavy ion collisions experiments performed at BNLs Relativistic Heavy Ion Collider 
(RHIC) create a hot and dense medium, the Quark-Gluon plasma (QGP). One of the most 
important discoveries at RHIC is the large elliptic flow in non-central Au+Au collisions, 
which is a clear indication of collective behavior. By now, it is generally agreed that the 
QGP thermalizes on times < 2.5 fm/c and behaves as a fluid with one of the lowest viscosity- 
to-entropy ratio observed in nature rj/s < 0.5 jl-3|. 

In recent years, relativistic hydrodynamics has become an efficient tool for describing the 
evolution of the fireball created at RHIC (there is a vast literature on this subject, see for 

nnn 

instance Refs. |4H15|). Ideal hydrodynamics has been partly successful in explaining the 



observed collective flow at low transverse momentum and in central collisions 



16| . Never- 



theless, one should notice that when using a realistic equation of state (including a crossover 
phase transition) and allowing for separate kinetic and chemical freeze outs, it seems difficult 
to fit the data with ideal hydrodynamics Moreover, if one aims eventually to de- 

rive the QGP viscosity from experimental data, one must start from a theoretical framework 
which allows for such effects. 

When one attempts to formulate a relativistic real hydrodynamics one finds that there 
is simply no equivalent to the non-relativistic Navier-Stokes equations. A straightforward 
relativistic generalization of the Navier-Stokes equations yields the so-called first order the- 

22I. One 



ories. These theories are plagued with causality and stability problems 
is therefore led to consider the so-called second-order theories (SOTs). These theories are 
presented as an expansion of the viscous tensor in velocity gradients, neglecting all orders 
higher than the second. They are unreliable in situations where these gradients are strong, 
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24]. 



and indeed they are known to fail, for example, in the description of strong shocks 
It is then valuable to develop alernative theories, not limited to weak velocity gradients, to 
provide at least an estimate of the expected accuracy of the gradient expansion. 



With this in mind, in Ref. 



25] the present authors developed an hydrodynamical de- 



scription of a conformal field within the framework of the so-called divergence-type theories 
(DTT) developed by Geroch [26[ (see also Refs. |27H29|). DTTs do not rely on velocity 



gradient expansions and in this sense they go beyond second-order theories. The purpose 
of this work is to present numerical results obtained from solving the equations of the DTT 
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in 2 + 1 dimensions. We use the equations to simulate Au+Au collisions, and compare the 
results both to experimental data and to a representative SOT. As was done in previous 
studies by other groups, we limit ourselves to boost-invariant longitudinal expansion in flat 
space-time. 

In the last years, there have been numerous theoretical studies of viscous hydrodynamics 
in 2 + 1 dimensions applied to heavy ion collisions 
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21 



30 



33| . Most previous works 



employ Israel- Stewart formalism or some variation of it in order to simulate the evolution of 

. I2I ] perfomed detailed simulations of 



the fireball. Very recently, Luzum and Romatschke [20|, 
Au+Au collisions based on the conformal hydrodynamical equations. The consistent picture 
that emerges from these diverse studies is that it is possible to match viscous hydrodynamics 
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2l|. See also [3J, |35| 



36] 



results to experimental data, provided rj/s < 5-6 x 1/47T 

In this paper we shall choose as prototype SOT the one developed by Baier et al 
(see also Bhattacharyya et al 37| and Natsuume et al 38j|). For simplicity we will refer to 
this hydrodynamic theory as BRSS (Baier- Romatschke-Son-Starinets). The BRSS theory 
is based on conformal invariance and extends the well-known Israel-Stewart (IS) formalism 
221 ] in that it contains all second-order terms that can appear in the stress-energy tensor of 



0, 



a conformal fluid [6j, [39H ■ The reason to study conformal field hydrodynamics is that it is 
relevant to the QGP since, as shown by Lattice calculations Q], QCD is approximately con- 
formal at high temperatures. In addition, a wealth of information of the strongly-interacting 
plasma such as transport coefficients (inaccesible to kinetic theory) can be obtained from 
the AdS/CFT correspondence 
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36 



38]. 



The main results we arrive at are: (i) the momentum anisotropy is smaller in the DTT 
than in the BRSS, (ii) the charged-hadron minumum-bias elliptic flow reaches its maximum 
value at lower pt in the DTT, and (iii) the matching of DTT results to data allows a 
viscosity-to-entropy ratio slightly larger than that of the BRSS. 

We should note here that in any hydrodynamic simulation there are numerous sources of 
uncertainty, such as those coming from the initial conditions, from the freeze-out procedure 
and from hadron dynamics (just to mention a few), which unfortunately prevent a precise 
determination of r]/s. Ours is not the exception, and here we merely intend to show that the 
DTT is an alternative to SOTs in the modeling of heavy-ion collisions (an alternative which 
may prove useful in those cases where large velocity gradients are present - e.g. shock- 
waves 23j). The results we obtain also show that the differences between hydrodynamic 
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formalisms are a significant source of uncertainty in the precise extraction of r]/s from data. 
This indubitably points to the conclusion that, until these uncertainties are under control, 
care should be taken when attempting to extract r]/s from hydrodynamic simulations. In 
this sense, the values for rj/s presented in this work should be regarded as rough estimates. 

The paper is organized as follows. In section |TT] we present the 2 + 1 hydrodynamic 
equations of the BRSS and of the DTT, and describe the initial conditions, the equation of 
state (EoS) and the freeze-out prescription employed in the simulations. In section II III we 
present and discuss the results obtained, and in section [TV] we present our conclusions. In 
Appendix |X] we evaluate the sensitivity of the results on the values of second-order transport 
coefficients and on the spatial mesh used in the simulations. In Appendix [B] we give a brief 
overview of divergence- type theories and derive Eq. (jSJ). 

II. THEORETICAL SETUP 
A. Hydrodynamic equations 

In this section we present the hydrodynamic equations of the BRSS and of the DTT for 
boost-invariant flow in 2 + 1 dimensions. We employ Milne coordinates defined by proper 
time t = \/t 2 — z 2 and rapidity ip = arctanh(2:/t), and, as mentioned in the Introduction, 
work in flat space-time. In these coordinates the metric tensor reads g^ v = (1, —1, — 1, — r 2 ). 
It is convenient to use Cartesian coordinates (x, y) in the transverse plane (instead of polar 
coordinates) since in this way the only non-vanishing Christoffel symbols are = r and 
rf^ = l/r. The fluid velocity is u = (u T , u x ,u y , 0) and is normalized as = 1. 

The stress-energy tensor for dissipative relativistic hydrodynamics is 

= prAt" - pA^ u + IP 1 ' with 

(1) 

A"" = g" u - v?u v 

where p and p are the energy density and the pressure in the local rest frame, and IP" is the 
viscous shear tensor which is transverse (w^EP" = 0), traceless and symmetric. The tensor 
A Mt/ is the spatial projector orthogonal to u** '. For a conformal fluid we have = 0, so 
p = 3p and the bulk viscosity vanishes. 

In what follows, Latin indices stand for transverse coordinates (x, y), is the geometric 
covariant derivative, D = u^D 11 and = A fl ' / D l , are the comoving time and space deriva- 
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tives, respectively, and < . . . > denote the spatial, symmetric and traceless projection of a 
tensor: 

A <„u> = (_ A ^a A ^u + 1 A ^ A av _ \ A ^ A ^A a ^ . (2) 

The hydrodynamic equations are the conservation equations for the stress-energy tensor 
together with the evolution equation for the shear tensor IP V . The former reads: 



(p + p)Du i = ^(g ij d jP - u l u a d aP ) - A^DpU^ 



D^-(/)+p)V/ + n"V 



where 



7/ ■ ?/ ■ 

_Dflii = ir a r — + — <9 r rr a + ^ir 



In the BRSS, the evolution of the shear tensor is given by [36|] 
d T U ia = — — rr a V u w M —IF* + — — — o 



(3) 



(4) 



3u r t 7T u t r 1T u T 

\— n< l n a> ^ £ ±Du» 5 

— &rr* 

where rj is the shear viscosity, (r n , Ai) are second-order transport coefficients, 

^ = v < ^m i/> (6) 

is the first-order shear tensor, and 

v T 

V H u" = <9 r w r + d^ 1 + — 

T 

V <x u x> = A™d T u x + A ix diU x - -A^Vpu" 

V <x u y> = ^A™d rU y + X -A^d T u x + ^A ix diU y (7) 

+ ^A iy d iU x - -A^VX 
V<^> = r 4 A^r^ - Va^VX . 
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We note that in Eq. (J5J) we have neglected terms involving the fluid vorticitv which would 



be multiplied by additional second-order transport coefficient s an d A 
is that, for two dimensional flow, it can be shown (see Ref. 



20 



36 



371 ] . The reason 



21] and references therein) 



that if the vorticity is zero initially (as it is the case here), it will remain negligible throughout 
the evolution up to terms which are third-order in velocity gradients (therefore beyond the 
scope of second-order theory). As already mentioned in the Introduction, the IS formalism 
is contained in the BRSS equations, as can be seen by setting Ai=0. 

In a DTT, the description of nonequilibrium hydrodynamic states requires the introduc- 



tion of a new tensor £ Q7 which is symmetric, traceless and vanishes in equilibrium 



26|. For 



a conformal fluid, C ai must be transverse as well. The DTT provides an equation of motion 



for ^ (see Ref. 
given by 



25] for details). For a conformal fluid in 2 + 1 dimensions the evolution is 



3r n f]u T 



The shear tensor is calculated from the nonequilibrium tensor £ Q7 as follows (see App. 

U ,u = _ h^El^a^ _ * A^C 7 6* 7 ) • (9) 

The transport coefficients of the BRSS and the DTT are the same because the DTT goes 
over to BRSS at second-order in velocity gradients. Note however that we are ignoring 
(possible) higher order corrections to the transport coefficients of the DTT. 

As independent variables we choose (p,u x ,u y ,U XX ,U xy ,U yy ) for the BRSS and 
(p,u x ,u y ,l; xx ,l; xy ,l; yy ) for the DTT. The r component of the velocity follows from nor- 
malization, u T = ^1 + v? x + Uy, while the other nontrivial components of (and of 
follow from the transversality and tracelessness conditions. 

In order to solve the hydrodynamic equations, we employ the method described in Ref. 
41 j (see also Ref. [20 1 ) . The set of six coupled differential equations is cast into a linear 



system for the time derivatives of the independent variables. This linear system is solved 
using a finite difference method which is first-order accurate in the temporal grid spacing 
and second-order accurate in the spatial grid spacing. To be more precise, the derivatives 



are replaced by: 



a,/ \ f(x,y,T + 6r) - f(x,y,r) 
drf{x } y } T) = and 

OT 

f(x + 5x, y, t) - f(x - 5x, y, r) 
25x 



(10) 



d x f(x,y,r) 



and similarly for the derivative in y. We have made nontrivial tests on the code such as, for 
example, use uniform initial data in the 2 + 1 numerical code to recover + 1 dimensional 
results (already presented in Ref. [251]). 



B. Initial conditions and transport coefficients 

Solution of the hydrodynamic equations requires initial conditions for the six independent 
variables. For the initial transverse velocity and shear tensor we use u x = u y = 0, which 
implies vanishing initial vorticit y, a nd (U xx ,U xy ,U yy ) = or (^ xx ,^ Ey ,^ yy ) = 0. It has been 



shown in several works 
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30 



311 ] that the evolution of the shear tensor IP 1 ' is quite 
insensitive to the initialization values, the difference being appreciable only at very early 
times. We have verified that the elliptic flow show very little sensitivity to the initialization 
of the shear tensor as well. In what follows, we take the initialization time to be tq = 1 
fm/c. 



42j , in which 



The initial energy density profile is calculated using a simple Glauber model 
for impact parameter b we have 

p(r , x,y,b) = C x oT A {x + ~, y)T A (x - - , y) (11) 

where a = 40 mb, C is a constant chosen such that p(r , 0,0,0) corresponds to a given 
initialization temperature To (via the equation of state), and T A is the nuclear thickness 
function given by 

/oo 
8 A (x,y,z) dz . (12) 
-oo 

The function 8 A is the Woods-Saxon density distribution for gold nuclei 

8 A (x, y, z) = — (13) 

1 + exp[(|x| - Rq)/x] 

with x = (x,y,z), Rq = 6.4 fm and \ = 0.54 fm. The parameter 5q is chosen such that 
j g? 3 x 5 A (x) = 197 as appropriate for Au nuclei. We note that in all calculations we use a 



7 



7.5 fm x 7.5 fm transverse plane (we have used a 13 fm x 13 fm transverse plane and found 
that there is no significant change in our results). 

Unless otherwise stated, we use values for the second-order transport coefficients corre- 
sponding to a strongly-coupled M = 4 Super- Yang Mills (SYM) plasma 
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r n = 2(2-ln2)4q and A x = , 
v ' sT 2ttT ' 



(14) 



where s is the entropy density. We will show in Appendix [X] that our results depend only 
weakly on the precise value of second-order transport coefficients. 



C. Equation of State 

The set of hydrodynamic equations must be closed with an equation of state (EoS). Since 
we are interested in computing elliptic flow of the produced particles we must use an EoS 
including hadronization. We employ the EoS by Laine and Schroder 43] which connects a 
high-order weak-coupling perturbative QCD calculation at high temperatures to a hadron 
resonance gas at low temperatures, via an analytic crossover (as suggested by Lattice QCD 



calculations 



18 



44JI). We notice that this EoS is the same as that used in Ref. 1201 1. so that 

n 

the comparisons we make are meaninful (see also Ref. [18|]). 



D. Freeze-out 



In order to compute the elliptic flow of produced particles, the freeze-out process must 
be simulated. To do so, we use the corresponding modules of the UVH2+1 code which 



is described in detail in Ref. 



20 



2l|. For completeness, we give here a brief overview of 



the Cooper-Frye freeze-out prescription 



45} implemented in UVH2+1. For the isothermal 



freeze-out we use here, the conversion from hydrodynamic to particle degrees of freedom 
takes place in a three-dimensional hypersurface. The spectrum for a single on-shell particle 
with momentum p 11 = (E, p) and degeneracy d is 

AN d 



E 



d 3 p (2tt) 



Vil dW f{x^jf) 



(15) 



S 



where <iS M is the normal vector on the hypersurface, and f is t 



function, customarily given by Grad's ansatz (see 
f(x*,p> t ) = f (x*,p> i ) + 6f(x»',p>') 
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20 



re non- equilibrium distribution 



21 



39 



46l | for details) 



/o(//) + / (/,p^ 



2T 2 (p + p) 



(16) 



cxp( 



T 



1 + 



2T> + p) 



The approximation in the second line holds when p » T, and it is used in our simulations. 



It has been shown in Ref. 



20 



21] that the systematic error of this approximation is very 



small at low tranverse momentum p? < 2.5 GeV, so we do not expect our results to have a 
significant error coming from this approximation. 

The UVH2+1 freeze-out module calculates the spectra for particle resonances with masses 
up to 2 GeV and then determines the spectra of stable particles including feed-down con- 
tributions (4?]] • For this last step it uses the AZHYDRO package j^sj. In this paper we will 
focus on the minimum-bias elliptic flow coefficient v 2 at central rapidity, which is given by 

fdbb v ( P T,b)v 2 (pT,b) 



v 2 (pi 



(17) 



fdbb vo(p T ,b) 

where the coeffients vq and v 2 are related to the particle spectra (including feed-down con- 
tributions) by 

= v (p T , b)[l + 2d 2 (p T , b) cos(20)] (18) 

with (f) = &rct&n(p y /p x ) and pt = (p x +Py)^ 2 - We note that kinetic and chemical freeze-out 
occur at the same temperature, which represents a simplification of the real process (see for 
instance Ref. 49j and references therein). 

We may obtain an internal consistency check of the hydrodynamical approximation by 
computing elliptic flow ffl8|) with the left hand side given by Eq. f Tl~5|) with the distribution 
function ( TT6|) . or else with only the equilibrium distribution function. The difference between 
these two results is the so-called nonequilibrium contribution to elliptic flow 5v 2 . 

It has been shown by Song and Heinz jsi| , by Chaudhuri [ill . [|3o| and by Dusling and 
Teaney [30] that beyond px ~2-3 GeV the nonequilibrium correction to the momentum 
distribution function becomes comparable to the equilibrium contribution, thus rendering 
Grad's ansatz unreliable. A way to estimate the value of px at which this happens is to 
compare 6v 2 to the total v 2 . We note that this way of determining when does Grad's 



ansatz become unreliable may miss corrections to the particle distribution function which 
are independent of azimuthal angle. However, it still provides an estimate which gives an 
idea of how large are the nonequilibrium corrections to the distribution function. 

In Figure [T] we show bv^jvi as a function of tranverse momentum, calculated with the 
BRSS and the DTT with rj/s = 0.08 and 6 = 7 fm. In this calculation, we set Tj = 333 MeV 
and Tf = 140 MeV in both models. 

0.8 I 1 1 1 1 1 1 1 1 — -. 

0.7 - 




0.5 1 1.5 2 2.5 3 3.5 4 4.5 

Pt [GeV] 



FIG. 1: (Color online) Nonequilibrium contribution to the total elliptic flow —bv^jvi calculated 
with the BRSS and DTT with 77/s = 0.08 and 6=7 fm. The initial and freeze-out temperatures 
are Tj = 333 MeV and Tf = 140 MeV in both models. The horizontal line indicates the (estimated) 
breakdown of Grad's ansatz. 

It is seen that, both in the BRSS and in the DTT, the nonequilibrium contribution to v 2 is 
significant even at low pt- In the BRSS 5v2 is always smaller that in the DTT, indicating that 
in the latter dissipative corrections to the nonequilibrium distribution function are larger. 
The difference between 8V2/V2 calculated with the DTT and the BRSS grows with pr, in 
particular for > 3 GeV where 8V2/V2 in the DTT starts growing faster. We can estimate 
the value of p? at which Grad's ansatz becomes unreliable as that at which —81)2/1)2 ~ 0.5, 
which is indicated by a horizontal line in the figure. In the DTT the value of px at which 
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Grad's ansatz becomes unreliable turns out to be px ~ 3.3 GeV instead of pr = 4 GeV as 
inferred in the BRSS. 



III. RESULTS 

In this section we go over to our main goal in this paper, namely to compare the results 
obtained with the DTT to RHIC data and to the BRSS. We start by comparing the entropy 
production and the spatial and momentum anisotropies in both hydrodynamic models, and 
then proceed to compare particle multiplicity, < pt > and elliptic flow to data in order to 
constrain the values of rj/s. 

Figure |2] shows the evolution of the total entropy production a in a collision with impact 
parameter 6 = 7 fm, for the DTT and the BRSS with rj/s = 0.08, an initial temperature 
Tj = 333 MeV and a freeze-out temperature Tf = 140 MeV. The total entropy production 
is given by the following expressions 



vbrss = / dxdy (d^S^ERss = / dxdy and 



a DTT = J dxdy {d p S^)\ DTT = J 



2r]T 

dxdy 



(19) 



2T 

We note that in order for the comparison to make sense we have used the same T and 77/s 
for the DTT and the BRSS. It is seen that the entropy production is significantly larger in 
the DTT. At early (late) times the DTT produces ~ 25 % (~ 40 %) more entropy than the 
BRSS. However, at very early times r — t$ < 1 fm/c the entropy production is the same in 
both models. 

We now go over to compare the spatial (e^) and momentum (e p ) anisotropies, defined as 

< y 2 — x 2 > < t xx — T yy > 

e x = n n — - an d Cv = — ^ — , (20) 

< y 2 + x 2 > p p <T XX + Tvy > y J 

where < ... > p denotes an average procedure over the tranverse plane with the energy 
density p as weighting factor. In Figure [3] we show e p calculated with the same value of 
v / s = 0.08, Ti =333 MeV and T f = 140 MeV for the BRSS and the DTT. This allows us 
to determine the differences between both hydrodynamic formalisms. It is seen that e p is 
slightly smaller in the DTT, indicating that the DTT leads to larger shear stress and thus to 
larger dissipation. This is consistent with the fact that, as already discussed, more entropy 
is produced in the DTT. 

11 
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FIG. 2: (Color online) Total entropy production a as a function of r — To for the BRSS and the 
DTT with r]js = 0.08, T t = 333 MeV, T f = 140 MeV and b = 7 fm. 

We now allow the values of rj/s to vary. In Figure H] we show the evolution of e x and 
e p for the BRSS and the DTT with b = 7 fm, T { =333 MeV, 7) = 140 MeV and several 
values of r]/s. The values of rj/s correspond to those that will be suitable for matching Kaon 
multiplicity and < px > to data, as will be shown later. 

It is seen that the spatial anisotropies calculated in the DTT and in the BRSS are very 
similar, and quite independent of rj/s. In contrast, the momentum anisotropies are seen to 
strongly depend on rj/s and are very different in the BRSS and the DTT. The DTT e p is 
systematically lower that the BRSS one, and, as we will see, this will be reflected on the 
calculated elliptic flow. 

In order to compare the elliptic flow calculated with the DTT or the BRSS to data and in 
this way constrain the value of r//s, we will follow the procedure of Luzum and Romatschke 
in j^. The idea is to determine the initial temperature Tj for b = and Tf by matching 
the hydrodynamic simulation to total multiplicity and < pt >■ As explained in Ref. (20 ]. 
one should refrain from trying to match pion multiplicity and < px > since the Boltzmann's 
approximation used in Eq. (|T6|) leads to an unavoidable systematic error. We will therefore 
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FIG. 3: (Color online) Evolution of the momentum anisotropy for the BRSS and the DTT with 
b = 7 fm, Ti = 333 MeV, Tj = 140 MeV and rj/s = 0.08. 

aim at a reasonable fit of Kaon multiplicity and < px > to determine Tj and Tf. In Appendix 
lAl we will show that the results for elliptic flow are only weakly dependent on the values of 
second-order transport coefficients (r n , Ai). For this reason, in the remaining of this section 
we will use coefficients appropiate for a SYM plasma as given in Eq. (j!4p . 

We start by presenting the results for total Kaon multiplicity and < pr > obtained in 
both models with rj/s = 0.08, T { = 333 MeV and T f = 140 MeV, compared to PHENIX 
data [51]. As in the case of spatial and momentum anisotropies, this allows us to evaluate 
the differences in both hydrodynamic formalisms for the same values of rj/s, Ti and Tf. We 
also show the results for the DTT with rj/s — 0.12 and T = 328 MeV (which as shown later 
is the value of T that gives, for this value of rj/s in the DTT, the best matching to data). 

Figure [5] shows the total multiplicity. The two sets of datapoints correspond to K + and 
K~ , which can not be distinguished in our model since the EoS corresponds to zero net- 
baryon density. It is seen that, in order to achieve the same multiplicity as in the BRSS, 
we must take rj/s = 0.12 in the DTT. Figure |6] shows < pr > obtained in both models. 
It is seen that < pt > is larger in the DTT than in the BRSS, indicating that the DTT 
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FIG. 4: (Color online) Evolution of the spatial (upper panel) and momentum (lower panel) 
anisotropies for the BRSS and the DTT with 6 = 7 fm, T = 333 MeV, T f = 140 MeV and 
different values of r\js. 

leads to larger transverse flow. Larger shear stress leads to larger transverse flow, so the 
results obtained for < px > confirm the fact that the DTT leads to stronger dissipation and 
consequently, larger entropy production. 

We now proceed to find Tj by matching total Kaon multiplicity and < pt > t o data. We 
note that we fix Tf = 140 MeV in both models. Figures [7] and |8] show the total multiplicity 
and < pt > calculated in the BRSS and the DTT. By performing this combined matching 
to multiplicity and < p T > data, we find for the DTT that T t = 328 MeV (323 MeV) at 
r]/s =0.12 (0.18) gives a reasonable fit to data, comparable to that obtained by the BRSS 



with % = 333 MeV (327 MeV) at rj/s =0.08 (0.16) |20|. 

With the values for Tj and Tf obtained before, we now go over to calculate the charged- 
hadron minimum-bias elliptic flow with the DTT and the BRSS and compare the results to 
data. In Figure M we show the elliptic flow calculated in the DTT and the BRSS with rj/s = 
0.08 and % = 333 MeV, and in the DTT with r]/s = 0.12 and % = 328 MeV, compared 



to experimental data from the STAR Collaboration 521 ] . We note that in order to estimate 
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FIG. 5: (Color online) Centrality dependence of total multiplicity for Kaons for Au+Au collisions 
at yfs = 200 GeV compared to the DTT and the BRSS for rj/s =0.08, T { = 333 MeV and T f = 



140 MeV. We also show the DTT resu 
is from the PHENIX Collaboration 



t with rj/s = 0.12 and T, = 328 MeV for comparison. Data 



5 it ] . The two sets of datapoints correspond to K + and K 



the removal 
(see Refs. 



: nonflow contributions to the elliptic flow, we reduce the STAR data by 20% 



20 



52]). It is seen that using the DTT with rj/s — 0.08 provides a poor fit to 



data, specially at intermediate and large transverse momentum. 

Figure [10] shows the charged-hadron minimum-bias elliptic of the BRSS and the DTT for 
different values of rj/s, compared to data. The most important differences in the calculated 
v 2 with the BRSS and with the DTT take place for p? > 2 GeV. It is seen that i>2 reaches 
its maximum at lower values of px in the DTT, and that the values for rj/s for which the 
DTT model is consistent with data {rj/s < 0.18) are slightly larger than those for the BRSS. 
The fact that the values of rj/s for which the elliptic flow is consistent with data are the 
same as those obtained from matching Kaon multiplicity and < pt > does not mean that 
t>2 measurements do not constrain rj/s, but rather that viscous hydrodynamics provides a 
consistent description of these observables in Au+Au collisions. 

We now wish to constrain the value of rj/s and give an upper bound beyond which our 
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FIG. 6: (Color online) Centrality dependence of < px > for Kaons for Au+Au collisions at y/s = 
200 GeV compared to the DTT and the BRSS for rj/s =0.08, T { = 333 MeV and T f = 140 MeV. 



We also show the DTT resu 
the PHENIX Collaboration 



t with rj/s = 0.12 and = 328 MeV for comparison. Data is from 



results cease to be consistent with experimental data. From Figures [5lfT0l we find that we 
can match the DTT results to experimental data provided rj/s < 0.18. It is difficult to 
determine the uncertainty in t]/s coming from the hydrodynamic simulation. Considering 
the dependence of final results on the following factors (i) the precise value of second-order 
transport coefficients (see Appendix [A]), (ii) the mesh grid (see Appendix [A} and (iii) the 
precise value of Tj and Tf obtained from matching hydrodynamic results to data on Kaon 
multiplicity and < pr >, as sources of uncertainty in the determination of rj/s, we can 
estimate this theoretical uncertainty in ± 0.07. Taking into account the uncertainty in the 
removal of non-flow contributions to the measured charged-hadron u 2 (min. bias) by STAR 



52| | . we can, following Ref. [20|], estimate the experimental uncertainty in ± 0.1. Therefore, 
we conclude that the DTT model favors rj/s < 0.35. We emphasize that our estimate for 
rj/s does not contemplate the uncertainty coming from several other factors such as bulk 



viscosity [53l. l55l. l56| . different temperatures for kinetic and chemical freeze-out [49|], precise 
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FIG. 7: (Color online) Centrality dependence of total multiplicity for Kaons for Au+Au collisions 
at y/s = 200 GeV compared to the DTT and the BRSS for various values of rj/s. The freeze-out 



temperature is Tf = 140 MeV. Data is from the PHENIX Collaboration 



5l| 



knowledge of the EoS and of the initial conditions [18|, |20|, l21|, [3lJ , which are expected to 
have a significant influence on the value of 7]/s. For this reason, our estimate for rj/s should 
be regarded as a conservative one. One should note, however, that both the DTT and the 
BRSS are consistent with experimental data, although none of the models can reproduce 
the saturation of the measured minimum bias elliptic flow. Recent studies |54j-l57j suggest 



that the origin of this failure of viscous hydrodynamics to reproduce saturation of the ellipic 
flow is Grad's quadratic ansatz for the nonequilibrium correction to the thermal distribution 
function given by Eq. (|T6|) . 



Finally, we note that t 
in several other works 



lis result for rj/s is in good agreement with the upper bound found 



20 



21| by similar matching of viscous hydrodynamics to data, and 



supports the notion that the matter created at RHIC exhibits almost perfect fluidity. 
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FIG. 8: (Color online) Centrality dependence of < px > for Kaons for Au+Au collisions at y/s = 
200 GeV compared to the DTT and the BRSS for various values of rj/s. The freeze-out temperature 



is T f = 140 MeV. Data is from PHENIX Collaboration 
IV. SUMMARY AND CONCLUSIONS 



We have studied the space-time evolution of a conformal plasma in 2 + 1 dimensions using 
second-order as well as divergence-type dissipative hydrodynamics. In the simulations, we 
employed a simple Glauber model to calculate the initial energy density distribution, a 
model equation of state with an analytic crossover, and the Cooper- Frye prescription for 



isothermal freeze-out (the latter implemented in the code UVH2+1 20|, l32|). 

We have made a comparison of the calculated Kaon total multiplicity and < pt > with 
experimental data by the PHENIX Collaboratio n |5ll| . as well as of the elliptic flow with 
experimental data by the STAR Collaboration [52]. We have found that the difference 
between the BRSS and the DTT elliptic flows starts to become significant when pt > 2 
GeV: the elliptic flow calculated with the DTT reaches its maximum at lower values of 
Pt- Including an estimate for the uncertainty in the determination of rj/s from data, we 
find that the DTT can be matched to RHIC data provided rj/s < 0.35, in good agreement 
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FIG. 9: (Color online) Comparison of STAR data on charged-hadron minimum-bias elliptic flow 
to DTT and BRSS results with tj/s = 0.08 and T, t = 333 MeV. We also show results of the DTT 
with rj/s = 0.12 and Tj = 328 MeV for comparison. 

with previous studies based on Israel-Stewart or BRSS equations. The results we obtain 
also show that the differences between hydrodynamic formalisms are a significant source of 
uncertainty in the precise extraction of T]/s from data. 

We note that niether the DTT nor the BRSS are able to reproduce the experimental 
saturation of elliptic flow at pt > 2.5 GeV, possibly pointing to the incorrectness of Grad's 
ansatz for the nonequilibrium distribution function (Eq. f|T6|) ). A related aspect that surely 
deserves further investigation is the inclusion of bulk viscosity in the hadronic stage of the 
fireball's evolution and during the deconfinement crossover. In this respect, it should be 



noted that recent work 55] has raised concern about the validity of Grad's ansatz for a 



reliable computation of freeze-out when bulk viscosity is present (see also Ref. 56[ for a 



different treatment of bulk viscosity). It was found in Ref. 55[ that the nonequilibrium 
corrections due to bulk viscosity to the distribution function are considerably larger than 
those coming from shear viscosity, rendering the application of Grad's moment method 
doubtful when bulk viscosity is taken into account. Moreover, Luzum and Ollitrault 
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521 ] . STAR event-plane data has been reduced 



FIG. 10: (Color online) Comparison of DTT resu 
minumum-bias elliptic flow by STAR (event-plane) 
by 20 % (red squares) to account for non-flow contributions (estimatively) jsj ]. We also show 
results for the BRSS for comparison. 

have found that even for a conformal plasma Grad's ansatz is disfavored by data on t> 4 / (V2) 2 , 
while Dusling, Moore and Teaney [53] have calculated the momentum dependence of the 
nonequilibrium contribution to the distribution function in a weak coupling setting and 
found it proportional to fi^ 2 . 

It would be therefore interesting to study this point within the framework of DTTs. In 
order to further investigate this issue, the relation between the DTT and microscopic theory 
must be precisely determined. 
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Appendix A: Dependence on transport coefficients and grid 

In this Appendix we will evaluate the dependence of the results obtained with the DTT on 
the spatial mesh and on the values for second-order transport coefficients. The parameters 
used for this evaluation are rj/s — 0.08, r = 1 fm/c, Tj = 333 MeV and Tf = 140 MeV in 
both cases. 

In Figure [TTJ we show the evolution of spatial anisotropy e x for a collision with b = 7 f m 
and a time step At = 0.002 fm/c, for different values of the space grid Ax = 0.1, 0.06 and 
0.04 fm. It is seen that the dependence of e x on Ax is very small and only significant at late 
times. There is practically no difference between the results obtained with Ax = 0.06 fm 
and 0.04 fm, and for this reason we employ Ax = 0.06 fm in all the simulations shown in 
this work. We have also checked that diminishing the time step did not change the results 
appreciably, and therefore use At = 0.002 fm/c throughout. 

0.35 | 1 1 1 1 1 1 1 1 1 1 




-0.05 1 1 1 1 1 1 1 1 1 1 1 

1 23456789 10 11 

T [fm/c] 

FIG. 11: (Color online) Evolution of the spatial anisotropy in the DTT with rj/s = 0.08 and b = 
7 fm, for different values of the space grid Ax. 
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A convenient quantity to measure the influence of second-order coefficients on final hadron 
observables is the elliptic flow coefficient v>2 (see Eq. (JT8"j) ). In Figure [12] we show the 
charged-hadron elliptic flow calculated in the DTT with rj/s — 0.08 and 6 = 7 fm, for two 
sets of second-order coefficients, namely those of a SYM plasma given in Eq. (TH|) and those 
corresponding to weakly-coupled Israel- Stewart formalism (r^ = 6r)/sT and Ai = 0). It is 
seen that the difference between v-i calculated with both sets of coefficients is negligable for 
Pt < 1.5 GeV and small for higher px- It is interesting to note that, although both curves 
are very similar, the elliptic flow calculated for the SYM plasma reaches its maximum at a 
lower value of px- We conclude that the results depend only weakly on the precise values of 
t w and Ai, at least for low values of rj/s. This is in agreement with the findings of Ref. [20 ]. 



0.15 




0.05 



FIG. 12: (Color online) Charged-hadron elliptic flow calculated in the DTT with rj/s = 0.08 and 
b = 7 fm with t w = 2(2 — ln2)r//sT and Ai = t]/2ttT, corresponding to a SYM plasma, and with 
TV = 6r//sT and Ai = 0, corresponding to weakly-coupled Israel-Stewart theory. 
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Appendix B: Divergence- type theories 



In this Appendix we give a brief summary of divergence-type theories (DTTs). Detailed 



discussions can be found in Refs. 



26 



28 



29). 



According to Geroch and Lindblom [26) , the hydrodynamical description of a nonequilib- 
rium state requires, besides the particle current and the stress-energy tensor T^, a new 
third order tensor A^ up obeying an equation of motion of divergence type. The dynamical 
equations are the conservation laws of and T^ V) together with an equation describing the 
dissipative part: 

D^ vp = I up (Bl) 

where A w and I up are algebraic local functions of iV M and T^ v and symmetric in the indices 
(z/, p), and is the covariant derivative. The entropy current is extended to 

S» = ^ - f3 u T^ u - aN» - A» vp i up (B2) 

where fi v = u v /T is the temperature vector, a = \xjT is the affinity, $ M is the thermodynamic 
potential and £ up is symmetric, traceless and vanish in equilibrium. 

We now require that the entropy and the thermodynamical potential be algebraic func- 
tions of (a, Pp, £iu,)- If the entropy production is to be nonnegative, then 

d&v <9$^ <9$^ 

-7T = N ^ ^r = TMl/ ; i^ = Aiivp ( B3 ) 

da dp v di V p 

Thus, as a consequence of the equations of motion, the entropy production rate is 

d^s* = -r p i vp . (B4) 

Since the stress-energy tensor is symmetric, we must also have 

*" = !; (B5) 

where is the so-called generating function of the theory. This means that 

every DTT is completely determined once x an d I are specified as algebraic functions of 
a, j3p, ^n V . The theory thus constructed satisfies the principles of relativity and entropy, and 



fully exploits the latter [28] 
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Introducing the symbol ( A to denote the set (a, A B the set (N^, T^, A flup ) and 
1b the set (0, 0, 7^), the theory is summed up in the equations 



/1M 

= -I B ( B 
Ib ■ 

We will now review the main results of Ref. 



(B6) 



D,A B 



25| where a quadratic DTT for a conformal 



fluid in d — 4 dimensions was developed. The most general generating function \ which is 
quadratic in ^ v is 



1 



X = X0(T) + X l{T)i^U V + —{Au^tLaTV? + B^pa 



(B7) 



where (a, A, £?) are coefficients to be determined. 

For a conformal field, T^ u — > e & ^( x ~ / )j 1 f lv under a Weyl transformation g^ v — >■ e~' 2w ( x ~ < >g llv , 
and g^T^ =0. These constraints pose no problem for Xo(^) an d XiCO> an d i n combination 
with the second equation of (IB3[) lead to 

a. 



Xo(T) 
Xi(T) 



6 
V 



-T 



T and 



(B8) 



where the energy density is p = aT 4 . 

For the quadratic part of x ? it is necessary to redefine the temperature in order to satisfy 
both constraints. To ensure that T^ v has the correct conformal weight we need a =6. The 
quadratic stress-energy tensor obtained from x then reads 

1 



rpflV 



-T" 4 [B{30u' i u u - 6AH£ pCT U + ZM^tp] ■ 



To get T% = we need A = 24B. The idea is to put 



rriflU 

l 2 



W v + 5Tlt v 



with W J,U traceless and transverse, and 5T% traceless. We get 



nf = AT' 



and 



(B9) 



(BIO) 



(Bll) 



(B12) 
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Defining the physical temperature T p by IT = T 4 + y — (5A/2a)£ pa £ p(7 we obtain 

T** = PpuV + ppA^ + (B13) 

with 

= r)^ v + AT' 4 - ^A pu ^^ (B14) 

and p p = aT 4 . 

The source term is written as I pu = + Iff", where If \ are linear and quadratic in 
£ pcr , respectively I± u is obtained by requiring the DTT to reproduce Eckart's theory at 
first order in velocity gradients, while I% v is obtained by requiring that the quadratic DTT 
satisfy the Second Law exactly The result is 

= + gT- 8 A^Z pu e a (B15) 

whereby the entropy production reads 

= ^Spotr • (B16) 

By requiring the DTT to reproduce the BRSS when is expanded at second-order in 
velocity gradients, the relation between (A,g) and (t t , Ai) is found to be 

A = — — L^T 8 (B17) 



and 



AlT ' (B18) 



9 

The equation of motion for the third-order tensor D p A pup = l vp renders an evolution 



equation for £ up '. In 2 + 1 the evolution is given by Eq. (jSJ). We note that in Ref. 25] we 
have set T = T p since the extra terms arising when using T p are of higher order (i.e. are 
terms which would be obtained from a cubic generating function x). We have numerically 
checked that these extra terms are negligible, thus one can use T = T p without appreciable 
change in the results. 
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